/*  eliptic_poly.c  */
/************************************************************************************
*																					*
*		Polynomial version of elliptic curve functions.  essentially the same as	*
*  optimal normal basis but uses polynomial functions.  Biggest difference is that	*
*  once poly_prime is picked (see polymain.c) you need to compute the "T matrix" to	*
*  help embed data onto a curve by solving a quadratic equation.					*
*																					*
*							Author = mike rosing									*
*							 date  = June 27, 1997									*
*																					*
************************************************************************************/

#include <stdio.h>
#include "field2n.h"
#include "poly.h"
#include "eliptic.h"

/*  Global parameters for polynomial math.  poly_prime is the irreducible polynomial
	for the Galois field arithmetic.  Tmatrix is used to solve quadratic equations
	and Smatrix is used to compute square roots.
*/

extern	FIELD2N poly_prime;
FIELD2N	Tmatrix[NUMBITS], Smatrix[NUMBITS], Trace_Vector;
INDEX	Null_Row;

/*  Subroutine to invert a binary matrix.  The method is brute force but simple.
	Start with a diagonal matrix in I, and for each row operation in mat_in that
	eliminates non diagonal terms, do the same in I.  The result is that mat_in 
	= 1 and I = 1/mat_in (transpose).  
	If the input transpose is set to any value, we then have to transpose I back
	into mat_out to get the correct result.  If transpose is zero, I is copied to
	mat_out and the value of error is returned to indicate any null row.
	
	Returns 0 if matrix invertible, row number of zero column if not invertible.
*/

INDEX poly_matrix_invert( mat_in, mat_out)
FIELD2N mat_in[NUMBITS], mat_out[NUMBITS];
{
	INDEX	row, col, i, j, rowdex, found, error;
	ELEMENT	src_mask, dst_mask;
	FIELD2N	dummy, Imatrix[NUMBITS];

/*  Create Imatrix as diagonalized  */

	for ( row=0; row<NUMBITS; row++)
	{
		null( &Imatrix[row]);
		Imatrix[row].e[NUMWORD - row/WORDSIZE] = 1L << (row%WORDSIZE);
	}
	error = 0;  /* hope this is return value  */
	
/*  Diagonalize input matrix.  Eliminate all other bits in each column. 
	First find a column bit that is set, and swap with diagonal
	row if needed.
*/

	for ( row = 0; row < NUMBITS; row ++)
	{
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);

/*  find a row of input matrix which has col = row bit set.
	First see if we get lucky, then do search.
*/
		found = 0;
		if ( !(mat_in[row].e[rowdex] & src_mask))
		{
			for (j = row+1; j<NUMBITS; j++)
			{
				if ( mat_in[j].e[rowdex] & src_mask)  /*  found one, swap rows  */
				{
					copy( &Imatrix[j], &dummy);
					copy( &Imatrix[row], &Imatrix[j]);
					copy( &dummy, &Imatrix[row]);
					copy( &mat_in[j], &dummy);
					copy( &mat_in[row], &mat_in[j]);
					copy( &dummy, &mat_in[row]);
					found = 1;
					break;
				}
			}
		}
		else found = 1;
		
/*  eliminate all other terms in this column  */

		if (found)
		{
			for ( i=0; i<NUMBITS; i++)
			{
				if ( i == row) continue;
				if ( mat_in[i].e[rowdex] & src_mask)
				{
					SUMLOOP(j)
					{
						Imatrix[i].e[j] ^= Imatrix[row].e[j];
						mat_in[i].e[j] ^= mat_in[row].e[j];
					}
				}
			}
		}			/*  end column eliminate  */
		else error = row;
	}				/*  end diagonalization   */


/*  Next step is to transpose diagonalized matrix.  This completes the 
	inversion process.  Clear Tmatrix to begin with. */
	
	for (row=0; row<NUMBITS; row++) null(&mat_out[row]);
	for (col = 0; col<NUMBITS; col++)
	{
		j = NUMWORD - col/WORDSIZE;
		src_mask = 1L << (col % WORDSIZE);
		for (row = 0; row<NUMBITS; row++)
		{
			if (Imatrix[row].e[j] & src_mask)
			{
				i = NUMWORD - row/WORDSIZE;
				dst_mask = 1L << (row  % WORDSIZE);
				mat_out[col].e[i] |= dst_mask;
			}
		}
	}
	return(error);
}

/*  initialize polynomial math for elliptic curve calculations.  
	
	Creates Tmatrix and Smatrix to help solve quadratic equations as
	well as Trace_Vector.  The Smatrix is only invoked to compute the
	square root of a polynomial modulo poly_prime.  The Tmatrix is a
	set of basis vectors which are summed assuming a solution to the
	quadratic exists.  The Trace_Vector is used to determine if the
	solution exists.  If an error occurs in creating the Smatrix the
	routine returns the row it could not diagonalize, otherwise the
	routine returns 0 for no errors.
*/

INDEX init_poly_math()
{
	INDEX	i, j, error, k, sum, row, rowdex, found, nulldex;
	ELEMENT	src_mask;
	FIELD2N	x, c, dummy;
	FIELD2N Trace[2*NUMBITS], Tz[NUMBITS], T2[NUMBITS];

/*  Create Trace_Vector for computing the trace in polynomial basis.  
	Given any input c = c_m*x^m + ... + c_1*x + c_0  the Trace_Vector
	will mask off and sum the correct coefficients of a_i.  The
	following method was suggested by Richard Pinch (rgep@dpmms.cam.ac.uk)
*/

/*  step 1 Build a list of powers of x modulo poly_prime from 0 to 2n-2 */

	null( &Trace[0]);
	Trace[0].e[NUMWORD] = 1L;
	for (i=1; i<2*NUMBITS-1; i++)
		mul_x_mod( &Trace[i-1], &poly_prime, &Trace[i]);

/*  step 2: Sum diagonals of nxn matrix using each row as a starting
	point.  Each sum amounts to the trace vector coefficient for that
	power of x.
*/

	null( &Trace_Vector);
	for (i=0; i<NUMBITS; i++)
	{
		sum = 0;
		j = NUMWORD;
		src_mask = 1;
		for ( k=i; k<i+NUMBITS; k++)
		{
			if ( Trace[k].e[j] & src_mask) sum ^= 1;
			src_mask <<= 1;
			if (!src_mask)
			{
				src_mask = 1;
				j--;
			}
		}
		if (sum) Trace_Vector.e[NUMWORD - i/WORDSIZE] |= 1L << (i % WORDSIZE);
	}
	
/*  Next compute Tz = z^2 + z matrix.  We already have every power of x in
	the Trace matrix, so each row of Tz is just two rows from Trace.  Do
	partial inversion to get Tmatrix.
*/
	for ( i=0; i<NUMBITS; i++)
	{
		SUMLOOP(j) 
		{
			Tz[i].e[j] = Trace[i].e[j] ^ Trace[i<<1].e[j];
			T2[i].e[j] = Trace[i<<1].e[j];
		}
	}
	
/*  Invert T2 matrix to get special square root matrix (called Smatrix)  */

	if (error = poly_matrix_invert( T2, Smatrix))
	{
		printf("Can not invert square root matrix.  Null row = %d\n", error);
		return(error);
	}

/*  Create a set of basis vectors for use in finding roots of z^2 + z = c
	This is similar to the matrix inversion, but there is no transpose and
	the row order is different.  Exactly how this works is not clear to me,
	but Prof. Pinch gets it and it works.
*/

/*  Create Tmatrix as diagonalized  */

	for ( row=0; row<NUMBITS; row++)
	{
		null( &Tmatrix[row]);
		Tmatrix[row].e[NUMWORD - row/WORDSIZE] = 1L << ( row % WORDSIZE);
	}

/*  Semi diagonalize input matrix.  
	First find a column bit that is set, and swap with diagonal
	row if needed.  Then eliminate all bits below it.
*/

	Null_Row = 0;
	nulldex = 0;
	for ( row = 0; row < NUMBITS; row++)
	{
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);
		
/*  find a row of input matrix which has col = row bit set.
	First see if we get lucky, then do search.
*/
		found = 0;
		if ( !(Tz[row].e[rowdex] & src_mask))
		{
			for (j = row+1; j<NUMBITS; j++)
			{
				if ( Tz[j].e[rowdex] & src_mask)  /*  found one, swap rows  */
				{
					copy( &Tmatrix[j], &dummy);
					copy( &Tmatrix[row], &Tmatrix[j]);
					copy( &dummy, &Tmatrix[row]);
					copy( &Tz[j], &dummy);
					copy( &Tz[row], &Tz[j]);
					copy( &dummy, &Tz[row]);
					found = 1;
			/*  keep track of original null row */
					if (row == nulldex) nulldex = j;
					break;
				}
			}
		}
		else found = 1;
		
/*  eliminate all other terms below diagonal in this column  */

		if (found)
		{
			for ( i=row+1; i<NUMBITS; i++)
			{
				if ( Tz[i].e[rowdex] & src_mask)
				{
					SUMLOOP(j)
					{
						Tmatrix[i].e[j] ^= Tmatrix[row].e[j];
						Tz[i].e[j] ^= Tz[row].e[j];
					}
				}
			}
		}			/*  end column eliminate  */
		else
		{		/*  mark null row and swap position with original  */
			Null_Row = row;
			copy( &Tmatrix[nulldex], &dummy);
			copy( &Tmatrix[row], &Tmatrix[nulldex]);
			copy( &dummy, &Tmatrix[row]);
			copy( &Tz[nulldex], &dummy);
			copy( &Tz[row], &Tz[nulldex]);
			copy( &dummy, &Tz[row]);
		}
	}				/*  end diagonalization   */

/*  finally eliminate all other terms above diagonal except for
	Null_Row.  Result is a set of basis vectors which converts
	c to z and z^2 + z = c.
*/

	for ( row = NUMBITS-1; row > 0; row--)
	{
		if (row == Null_Row) continue;
		
		rowdex = NUMWORD - row/WORDSIZE;
		src_mask = 1L << (row % WORDSIZE);

		for (i = row-1; i>=0; i--)
		{
			if (Tz[i].e[rowdex] & src_mask)
			{
				SUMLOOP(j)
				{
					Tmatrix[i].e[j] ^= Tmatrix[row].e[j];
					Tz[i].e[j] ^= Tz[row].e[j];
				}
			}
		}
	}			/*  end column eliminate  */
	return(0);
}

/*********************************************************************
*																	 *
*	solve a quadratic equation over an irreducible polynomial field. *
*   Enter with parameters for the equation y^2 + xy + f = 0, where y *
*  is the unknown, x is	usually the x coordinate of a point and 	 *
*  f = f(x) of a curve.  											 *
*																	 *
*	Assumes init_poly_math already run with Trace_Vector, Null_Row   *
*   and Tmatrix filled out.											 *
*																	 *
*	computes c = f/x^2 and tests if the trace condition is met.  If  *
*   not	returns an error code of 1.  Otherwise it means we can solve *
*   the reduced	equation z^2 + z + c = 0.  The change of variable is *
*   y = xz.															 *
*																	 *
*	returns error code zero, y[0] = xz, y[1] = y[0] + x				 *
*********************************************************************/

INDEX poly_quadratic( x, f, y)
FIELD2N *x, *f, y[2];
{
	FIELD2N	c, z, dummy;
	FIELD2N test1, test2;
	INDEX	i, j, k;
	ELEMENT	sum, mask;

/*  first check to see if x is zero  */

	sum = 0;
	SUMLOOP(i) sum |= x->e[i];
	
	if (sum)
	{
/*  compute c = x^-2 * f  */

		copy (x, &c);
		poly_mul(&c, x, &dummy);		/*  get x^2  */
		poly_inv(&dummy, &z);
		poly_mul( f, &z, &c);

/*  verify that trace of c = 0.  If not, no solution
	is possible.  */
	
		sum = 0;
		SUMLOOP(i) sum ^= c.e[i] & Trace_Vector.e[i];
        mask = ~0;
        for (i = WORDSIZE/2; i > 0; i >>= 1)
	    {
    	    mask >>= i;
        	sum = ((sum & mask) ^ (sum >> i));
        }
 
   	/* if last bit is set, there is no solution to equation.  
	This eliminates half the points in the field, which might
	make sense to a mathematician.
	*/
		if (sum)
		{
			null( &y[0]);
			null( &y[1]);
			return(1);
		}

/*  clear out null row bit.  that part of matrix will not work.  */

		j = NUMWORD - Null_Row / WORDSIZE;
		c.e[j] &= ~( 1L << (Null_Row % WORDSIZE));

/*  for every bit set in c, add that row of Tmatrix to solution.  */

		null( &z);
		mask = 1;
		j = NUMWORD;
		for (i=0; i<NUMBITS; i++)
		{
			if ( c.e[j] & mask)
				SUMLOOP(k) z.e[k] ^= Tmatrix[i].e[k];
			mask <<= 1;
			if ( !mask)
			{
				mask = 1;
				j--;
			}
		}
		
/*  compute final solution using input parameters.  */

		poly_mul( x, &z, &y[0]);	
		SUMLOOP(i) y[1].e[i] = y[0].e[i] ^ x->e[i];
		return(0);
	}
	else
	{
/*  x input was zero.  Return y = square root of f. Process
	involves ANDing each row of Smatrix with f and summing
	all bits to find coefficient for that power of x */

		null( &z);
		for (j = 0; j<NUMBITS; j++)
		{
			sum = 0;
			SUMLOOP(i) sum ^= Smatrix[j].e[i] & f->e[i];
			if (sum)
			{
        		mask = -1L;
        		for (i = WORDSIZE/2; i > 0; i >>= 1)
	       	 	{
    	    		mask >>= i;
        	    	sum = ((sum & mask) ^ (sum >> i));
            	}
	        } 
	        if (sum) z.e[NUMWORD - j/WORDSIZE] |= (1L << j%WORDSIZE);
		}
		copy( &z, &y[0]);
		copy( &z, &y[1]);
		return(0);
	}
}

/*  print field matrix.  an array of FIELD2N of length NUMBITS is assumed.  The
	bits are printed out as a 2D matrix with an extra space every 5 characters.
*/

void matrix_print( name, matrix, file)
FIELD2N matrix[NUMBITS];
char *name;
FILE *file;
{
	INDEX i,j;
	
	fprintf(file,"%s\n",name);
	for (i = NUMBITS-1; i>=0; i--)
	{
		if (i%5==0) fprintf(file,"\n");	/* extra line every 5 rows  */
		for (j = NUMBITS-1; j>=0; j--)
		{
			if ( matrix[i].e[NUMWORD - j/WORDSIZE] & (1L<<(j%WORDSIZE)))
				fprintf(file,"1");
			else fprintf(file,"0");
			if (j%5==0) fprintf(file," ");  /* extra space every 5 characters  */
		}
		fprintf(file,"\n"); 
	}
	fprintf(file,"\n");
}

/*  compute f(x) = x^3 + a_2*x^2 + a_6 for non-supersingular elliptic curves */

void poly_fofx( x, curv, f)
FIELD2N *x, *f;
CURVE *curv;
{
	FIELD2N 	x2, x3;
	INDEX		i;
	
	copy ( x, &x3);
	poly_mul( x, &x3, &x2);       /*  get x^2  */
	if (curv->form) poly_mul ( &x2, &curv->a2, f);
	else null( f);
	poly_mul( x, &x2, &x3);       /*  get x^3  */
	SUMLOOP (i) f->e[i] ^= ( x3.e[i] ^ curv->a6.e[i] );
}

/*  embed data onto a curve.
	Enter with data, curve, ELEMENT offset to be used as increment, and
	which root (0 or 1).
	Returns with point having data as x and correct y value for curve.
	Will use y[0] for last bit of root clear, y[1] for last bit of root set.
	if ELEMENT offset is out of range, default is 0.
*/

void poly_embed( data, curv, incrmt, root, pnt)
FIELD2N	*data;
CURVE	*curv;
INDEX	incrmt, root;
POINT	*pnt;
{
	FIELD2N		f, y[2];
	INDEX		inc = incrmt;
	INDEX		i;
	
	if ( (inc < 0) || (inc > NUMWORD) ) inc = 0;
	copy( data, &pnt->x);
	poly_fofx( &pnt->x, curv, &f);
	while (poly_quadratic( &pnt->x, &f, y))
	{
		pnt->x.e[inc]++;
		poly_fofx( &pnt->x, curv, &f);
	}
	copy ( &y[root&1], &pnt->y);
}

/****************************************************************************
*                                                                           *
*   Implement elliptic curve point addition for polynomial basis form.  	*
*  This follows R. Schroeppel, H. Orman, S. O'Mally, "Fast Key Exchange with*
*  Elliptic Curve Systems", CRYPTO '95, TR-95-03, Univ. of Arizona, Comp.   *
*  Science Dept.                                                            *
****************************************************************************/

void poly_esum (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    INDEX   i;
    FIELD2N  x1, y1, theta, onex, theta2;
    ELEMENT  check;
	
/*  check if p1 or p2 is point at infinity  */

	check = 0;
	SUMLOOP(i) check |= p1->x.e[i] | p1->y.e[i];
	if (!check)
	{
		printf("summing point at infinity\n");
		copy_point( p2, p3);
		return;
	}
	check = 0;
	SUMLOOP(i) check |= p2->x.e[i] | p2->y.e[i];
	if (!check) 
	{
		printf("summing point at infinity\n");
		copy_point( p1, p3);
		return;
	}
	
/*  compute theta = (y_1 + y_2)/(x_1 + x_2)  */

    null(&x1);
    null(&y1);
    check = 0;
    SUMLOOP(i) 
    {
		x1.e[i] = p1->x.e[i] ^ p2->x.e[i];
		y1.e[i] = p1->y.e[i] ^ p2->y.e[i];
		check |= x1.e[i];
    }
    if (!check)   /* return point at infinity  */
    {
    	printf("input to elliptic sum has duplicate x values\n"); 
    	null(&p3->x);
    	null(&p3->y);
    	return;
    }
    poly_inv( &x1, &onex);
    poly_mul( &onex, &y1, &theta);  /*  compute y_1/x_1 = theta */
    poly_mul(&theta, &theta, &theta2);   /* then theta^2  */

/*  with theta and theta^2, compute x_3  */

    if (curv->form)
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i];

/*  next find y_3  */

    SUMLOOP (i) x1.e[i] = p1->x.e[i] ^ p3->x.e[i];
    poly_mul( &x1, &theta, &theta2);
    SUMLOOP (i) p3->y.e[i] = theta2.e[i] ^ p3->x.e[i] ^ p1->y.e[i];
}

/*  elliptic curve doubling routine for Schroeppel's algorithm over polymomial
    basis.  Enter with p1, p3 as source and destination as well as curv
    to operate on.  Returns p3 = 2*p1.
*/

void poly_edbl (p1, p3, curv)
POINT *p1, *p3;
CURVE *curv;
{
    FIELD2N  x1, y1, theta, theta2, t1;
    INDEX   i;
    ELEMENT  check;
    
    check = 0;
    SUMLOOP (i) check |= p1->x.e[i];
    if (!check)
    {
    	printf("doubling point at infinity\n");
    	null(&p3->x);
    	null(&p3->y);
    	return;
    }

/*  first compute theta = x + y/x  */

    poly_inv( &p1->x, &x1);
    poly_mul( &x1, &p1->y, &y1);
    SUMLOOP (i) theta.e[i] = p1->x.e[i] ^ y1.e[i];

/*  next compute x_3  */

    poly_mul( &theta, &theta, &theta2);
    if(curv->form)
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i];

/*  and lastly y_3  */

    theta.e[NUMWORD] ^= 1;			/*  theta + 1 */
    poly_mul( &theta, &p3->x, &t1);
    poly_mul( &p1->x, &p1->x, &x1);
    SUMLOOP (i) p3->y.e[i] = x1.e[i] ^ t1.e[i];
}

/*  subtract two points on a curve.  just negates p2 and does a sum.
    Returns p3 = p1 - p2 over curv.
*/

void poly_esub (p1, p2, p3, curv)
POINT   *p1, *p2, *p3;
CURVE   *curv;
{
    POINT   negp;
    INDEX   i;

    copy ( &p2->x, &negp.x);
    null (&negp.y);
    SUMLOOP(i) negp.y.e[i] = p2->x.e[i] ^ p2->y.e[i];
    poly_esum (p1, &negp, p3, curv);
}

/*  need to move points around, not just values.  Optimize later.  */

void copy_point (p1, p2)
POINT *p1, *p2;
{
	copy (&p1->x, &p2->x);
	copy (&p1->y, &p2->y);
}

/*  Routine to compute kP where k is an integer (base 2, not normal basis)
	and P is a point on an elliptic curve.  This routine assumes that K
	is representable in the same bit field as x, y or z values of P.  Since
	the field size determines the largest possible order, this makes sense.
    Enter with: integer k, source point P, curve to compute over (curv) 
    Returns with: result point R.

  Reference: Koblitz, "CM-Curves with good Cryptografic Properties", 
	Springer-Verlag LNCS #576, p279 (pg 284 really), 1992
*/

void  poly_elptic_mul(k, p, r, curv)
FIELD2N	*k;
POINT	*p, *r;
CURVE	*curv;
{
	char		blncd[NUMBITS+1];
	INDEX		bit_count, i;
	ELEMENT		notzero;
	FIELD2N		number;
	POINT		temp;

/*  make sure input multiplier k is not zero.
	Return point at infinity if it is.
*/
	copy( k, &number);
	notzero = 0;
	SUMLOOP (i) notzero |= number.e[i];
	if (!notzero)
	{
		null (&r->x);
		null (&r->y);
		return;
	}

/*  convert integer k (number) to balanced representation.
	Called non-adjacent form in "An Improved Algorithm for
	Arithmetic on a Family of Elliptic Curves", J. Solinas
	CRYPTO '97. This follows algorithm 2 in that paper.
*/
	bit_count = 0;
	while (notzero)
	{
	/*  if number odd, create 1 or -1 from last 2 bits  */
	
		if ( number.e[NUMWORD] & 1 )
		{
			blncd[bit_count] = 2 - (number.e[NUMWORD] & 3);
			
	/*  if -1, then add 1 and propagate carry if needed  */
			
			if ( blncd[bit_count] < 0 )
			{
				for (i=NUMWORD; i>=0; i--)
				{
					number.e[i]++;
					if (number.e[i]) break;
				}
			}
		}
		else
			blncd[bit_count] = 0;
	
	/*  divide number by 2, increment bit counter, and see if done  */
	
		number.e[NUMWORD] &= ~0 << 1;
		rot_right( &number);
		bit_count++;
		notzero = 0;
		SUMLOOP (i) notzero |= number.e[i];
	}
		
/*  now follow balanced representation and compute kP  */

	bit_count--;
	copy_point(p,r);		/* first bit always set */
	while (bit_count > 0) 
	{
	  poly_edbl(r, &temp, curv);
	  bit_count--;
	  switch (blncd[bit_count]) 
	  {
	     case 1: poly_esum (p, &temp, r, curv);
				 break;
	     case -1: poly_esub (&temp, p, r, curv);
				  break;
	     case 0: copy_point (&temp, r);
	   }
	}
}

void print_field( string, x)
char *string;
FIELD2N *x;
{
	INDEX i;
	
	printf("%s\n",string);
	SUMLOOP(i) printf("%8x ",x->e[i]);
	printf("\n");
}

void print_point( string, point)
char *string;
POINT *point;
{
	INDEX i;
	
	printf("%s\n",string);
	printf("x: ");
	SUMLOOP(i) printf("%8x ",point->x.e[i]);
	printf("\n");
	printf("y: ");
	SUMLOOP(i) printf("%8x ",point->y.e[i]);
	printf("\n");
}

void print_curve( string, curv)
char *string;
CURVE *curv;
{
	INDEX i;
	
	printf("%s\n",string);
	printf("form: %d\n",curv->form);
	if (curv->form)
	{
		printf("a2: ");
		SUMLOOP(i) printf("%lx ",curv->a2.e[i]);
		printf("\n");
	}
	printf("a6: ");
	SUMLOOP(i) printf("%lx ",curv->a6.e[i]);
	printf("\n\n");
}

